
####---------------------------------------------------------#####
# This script reproduces all results reported in the Letter and  #
# in the Supplementary Materials. The matching and calculation   #
# of expected durations from the Cox models take a while to run. #
# The results from these procedures are therefore provided in    #
# R images, so that the code can be ran with these procedures    #
# commented out. To verify the matching and expected durations,  #
# please uncomment the relevant code.                            #
##################################################################

#### Loading  packages ####

library(MatchIt)
library(survival)
library(coxed)
library(ggplot2)
library(cobalt)
library(texreg)
library(ggplot2)
library(survminer)

library(simPH)
library(dplyr)

## working directory (replace path with your own directory) ###
setwd("C:/Users/oyvinsti/Dropbox/!PhD/PaperProjects/ISA2017/BJPSRevision/ReplicationCode") 
### Loading data ###
load("DirectingComplianceReplicationDataApril2019.RData")



## setting seed #####
set.seed(232321)

### Models estimated before matching:
model1 <- coxph(Surv(tstart, tstop, compliance) ~ 
                  article46tvc+
                  cluster(respondent),
                data = implementation)
summary(model1)
cox.zph(model1, "rank")


## this takes a while to run, commented out to save time:
# coxed1 <- coxed(model1,
#                 id = implementation$id,
#                 method = "npsf",
#                 newdata = mutate(implementation, article46tvc = 0),
#                 newdata2 = mutate(implementation, article46tvc = 1),
#                 bootstrap = TRUE,
#                 B = 1000)
# save(coxed1,
#      file = "coxed1.RData")
load("coxed1.RData")

summary(coxed1)




#### Selecting the first row for each case to use for the matching ####
firstRow <- implementation[which(implementation$tstart == 0),]


### Matching #####
matchingData <- na.omit(firstRow[,c('id',
                                    'article46judgmentAnyPoint',
                                    'legislation',
                                    'jurisprudence',
                                    'exe_adm',
                                    'publication',
                                    'practical',
                                    'prosecution',
                                    'property',
                                    'reopen',
                                    'other_individual', 
                                    'numberOfArticlesViolated',
                                    "article2violated",
                                    "article3violated",
                                    "article5violated",
                                    "article6violated",
                                    "article8violated",
                                    "article10violated",
                                    "article13violated",
                                    "article14violated",
                                    "protocol1_1violated",
                                    "Judgment.Year",
                                    "postProtocol11",
                                    "postCoMWorkingmethodChange",
                                    "postProtocol14",
                                    'nondemocracy',
                                    'v2x_accountability',
                                    'execap',
                                    'polconiii',
                                    'newdemocracy1')])





### this takes a while to run, commented out to save time:

# genetic.any <- matchit(formula = article46judgmentAnyPoint ~ legislation + jurisprudence +
#                          exe_adm + publication + practical + prosecution + property +
#                          reopen + other_individual + numberOfArticlesViolated + article2violated +
#                          article3violated + article5violated + article6violated +
#                          article8violated + article10violated + article13violated +
#                          article14violated + protocol1_1violated + polconiii + nondemocracy +
#                          newdemocracy1 + execap + v2x_accountability + Judgment.Year +
#                          postProtocol11 + postCoMWorkingmethodChange + postProtocol14,
#                        data = matchingData,
#                        method = "genetic",
#                        unif.seed = 232321,
#                        wait.generations = 10,
#                        nboots = 1000,
#                        ratio = 1,
#                        pop.size = 1000,
#                        replace = TRUE,
#                        ties = TRUE)
# save(genetic.any,
#      file = "firstMatching.RData")

load("firstMatching.RData") ## loading the results from the genetic matching

summary(genetic.any)




pdf("figure1.pdf")
love.plot(bal.tab(genetic.any), 
          var.names = data.frame(old =  c(
            'legislation',
            'jurisprudence',
            'exe_adm',
            'publication',
            'practical',
            'prosecution',
            'property',
            'reopen',
            'other_individual', 
            'numberOfArticlesViolated',
            "article2violated",
            "article3violated",
            "article5violated",
            "article6violated",
            "article8violated",
            "article10violated",
            "article13violated",
            "article14violated",
            "protocol1_1violated",
            "Judgment.Year",
            "postProtocol11",
            "postCoMWorkingmethodChange",
            "postProtocol14",
            'nondemocracy',
            'execap',
            'polconiii',
            'newdemocracy1',
            'v2x_accountability'), 
            new =  c(
              'Legislative change',
              'Jurisprudential change',
              'Executive action',
              'Publication and dissemination',
              'Practical measure',
              'Domestic investigation or prosecution',
              'Return of property',
              'Reopening of domestic proceedings',
              'Other individial measure', 
              'Number of articles violated',
              "Right to life violation",
              "Prohibition of torture violation",
              "Right to liberty violation",
              "Right to fair trial violation",
              "Right to privacy and family life violated",
              "Freedom of expression violation",
              "Right to effective remedy violation",
              "Prohibition of discrimination violation",
              "Property rights violations",
              "Judgment year",
              "After protocol 11",
              "After change in CoM working methods",
              "After protocol 14",
              'Non-democracy',
              'Bureaucratic capacity',
              'Political constraints',
              'New democracy',
              'Government accountability')), 
          drop.distance = TRUE,
          var.order = "unadjusted", 
          abs = FALSE, 
          stat = "mean.diffs",
          colors = c("black","grey"), 
          shapes = c(17,19))+
  theme_classic()
dev.off()




matchingData$genetic.any.weights  <- genetic.any$weights



implementation <- merge(implementation,
                        matchingData[,c("id","genetic.any.weights")],
                        by = "id",
                        all.x = TRUE)

model2  <- coxph(Surv(tstart, tstop, compliance) ~ 
                   article46tvc  +
                   +cluster(respondent),
                 weights = implementation[which(implementation$genetic.any.weights > 0),]$genetic.any.weights,
                 data = implementation[which(implementation$genetic.any.weights > 0),])
summary(model2)
cox.zph(model2, "rank")


model2Data <- implementation[which(implementation$genetic.any.weights > 0),]
# coxed2 <- coxed(model2,
#                 id = model2Data$id,
#                 method = "npsf",
#                 newdata = mutate(model2Data, article46tvc = 0),
#                 newdata2 = mutate(model2Data, article46tvc = 1),
#                 bootstrap = TRUE,
#                 B = 1000)
# save(coxed2,
#      file = "coxed2.RData")
# 

load("coxed2.RData")

summary(coxed2)



model3 <- coxph(Surv(tstart, tstop, compliance) ~ 
                  article46tvc+
                  legislation +
                  jurisprudence +
                  exe_adm +
                  publication +
                  practical +
                  property +
                  reopen +
                  other_individual +
                  numberOfArticlesViolated+
                  years_since_legislative_election_tvc +
                  followCases+
                  polconiii+
                  nondemocracy+
                  newdemocracy1+
                  execap +
                  v2x_accountability+
                  cluster(respondent),
                weights = implementation[which(implementation$genetic.any.weights > 0),]$genetic.any.weights,
                data = implementation[which(implementation$genetic.any.weights > 0),],
                x = TRUE)

cox.zph(model3, "rank")

summary(model3)

model3Data <- na.omit(model2Data[,c("tstart", "tstop", "compliance", "id", 
                                    "genetic.any.weights",
                                    "respondent",
                                    names(coefficients(model3)))])

# coxed3 <- coxed(model3,
#                 id = model3Data$id,
#                 method = "npsf",
#                 newdata = mutate(model3Data, article46tvc = 0),
#                 newdata2 = mutate(model3Data, article46tvc = 1),
#                 bootstrap = TRUE,
#                 B = 1000)
# save(coxed3,
#      file = "coxed3.RData")
load("coxed3.RData")

summary(coxed3)

model4Data <- model3Data
model4Data$article46tvcXv2x_accountability <- model4Data$article46tvc *  model4Data$v2x_accountability
model4 <- coxph(Surv(tstart, tstop, compliance) ~ 
                  article46tvc +  v2x_accountability +
                  article46tvcXv2x_accountability +
                  legislation +
                  jurisprudence +
                  exe_adm +
                  publication +
                  practical +
                  property +
                  reopen +
                  other_individual +
                  numberOfArticlesViolated+
                  years_since_legislative_election_tvc +
                  followCases+
                  polconiii+
                  nondemocracy+
                  newdemocracy1+
                  execap +
                  cluster(respondent),
                weights = model4Data$genetic.any.weights,
                data = model4Data,
                tt = logtime,
                x = TRUE)
cox.zph(model4, "rank")
summary(model4)



# coxed4high <- coxed(model4,
#                     id = model4Data$id,
#                     method = "npsf",
#                     newdata = mutate(model4Data, article46tvc = 0,
#                                      v2x_accountability = max(model4Data$v2x_accountability),
#                                      article46tvcXv2x_accountability = 0 * max(model4Data$v2x_accountability)),
#                     newdata2 = mutate(model4Data, article46tvc = 1,
#                                       v2x_accountability = max(model4Data$v2x_accountability),
#                                       article46tvcXv2x_accountability = 1 * max(model4Data$v2x_accountability)),
#                     bootstrap = TRUE,
#                     B = 1000)
# 
# save(coxed4high,
#      file = "coxed4high.RData")

load("coxed4high.RData")

summary(coxed4high)

# coxed4low <- coxed(model4,
#                    id = model4Data$id,
#                    method = "npsf",
#                    newdata = mutate(model4Data, article46tvc = 0,
#                                     v2x_accountability = min(model4Data$v2x_accountability),
#                                     article46tvcXv2x_accountability = 0 * min(model4Data$v2x_accountability)),
#                    newdata2 = mutate(model4Data, article46tvc = 1,
#                                      v2x_accountability = min(model4Data$v2x_accountability),
#                                      article46tvcXv2x_accountability = 1 * min(model4Data$v2x_accountability)),
#                    bootstrap = TRUE,
#                    B = 1000)
# save(coxed4low,
#      file = "coxed4low.RData")

load("coxed4low.RData")




# coxed4median <- coxed(model4,
#                       id = model4Data$id,
#                       method = "npsf",
#                       newdata = mutate(model4Data, article46tvc = 0,
#                                        v2x_accountability = median(model4Data$v2x_accountability),
#                                        article46tvcXv2x_accountability = 0 * median(model4Data$v2x_accountability)),
#                       newdata2 = mutate(model4Data, article46tvc = 1,
#                                         v2x_accountability = median(model4Data$v2x_accountability),
#                                         article46tvcXv2x_accountability = 1 * median(model4Data$v2x_accountability)),
#                       bootstrap = TRUE,
#                       B = 1000)
# save(coxed4median,
#      file = "coxed4median.RData")
load("coxed4median.RData")

summary(coxed4median)


coxeds <- rbind(coxed4high$median.diff,
                coxed4median$median.diff,
                coxed4low$median.diff,
                coxed3$median.diff,
                coxed2$median.diff,
                coxed1$median.diff )

coxeds <- coxeds / 365.25

coxeds$description <- c(   "6. Remedial indications, 
                        multivariate interaction model on matched dataset, 
                           highest level of government accountability",
                           "5. Remedial indications, 
                        multivariate interaction model on matched dataset, 
                           median level of government accountability",
                           "4. Remedial indications, 
                        multivariate interaction model on matched dataset, 
                           lowest level of government accountability",
                           "3. Remedial indications, 
                        multivariate model on matched dataset",
                           "2. Remedial indications, 
                        matched dataset",
                           "1. Remedial indications,
                        unadjusted")

coxeds$order <- 1:6

coxeds$model <- factor(coxeds$order, 
                       labels = coxeds$description)






pdf("figure2.pdf")
ggplot(coxeds)+
  geom_segment(aes(x = model, 
                   xend = model, 
                   y = lb, 
                   yend = ub))+
  geom_hline(aes(yintercept = 0), colour = "grey")+
  geom_point(aes(x = model, y = median))+
  coord_flip()+
  ylab(expression(paste(Delta, " Expected years until compliance")))+
  xlab("")+
  theme_classic()
dev.off()



### Regression table ####


mainModels <- list(model1, model2, model3, model4)
for(i in 1:length(mainModels)){
  mainModels[[i]]$n <- sum( mainModels[[i]]$y[,1] == 0)
}

texreg(mainModels , 
       fontsize = "scriptsize",
       label = "tab:FullCoxModels",
       caption = "Full Cox models", 
       file = "MainModels.tex",
       include.rsquared = FALSE,
       include.maxrs= FALSE, 
       caption.above = TRUE, 
       include.missings = FALSE,
       include.zph = FALSE,
       custom.coef.map = list(article46tvc = "Remedial indications",
                              v2x_accountability  = "Government accountability", 
                              article46tvcXv2x_accountability  = "Government accountability * Remedial indications", 
                              years_since_legislative_election_tvc  = "Years since last election", 
                              polconiii =  "Political constraints",
                              nondemocracy = "Non-democracy",
                              newdemocracy1 = "New democracy",
                              execap = "Bureaucratic capacity",
                              legislation =  "Need for legislative change",
                              jurisprudence = "Need for jurisprudential change",
                              exe_adm = "Need for executive action",
                              publication = "Need for publication/dissemination",
                              practical = "Need for practical measures",
                              property =  "Need for property return",
                              reopen =  "Need for reopening of domestic case",
                              other_individual =  "Need for other individual measure",
                              numberOfArticlesViolated =  "Number of articles violated",
                              followCases = "Repetitive cases"))



#### interaction plot for appendix #####

model4AlternativeCall <- model4$call$formula
model4AlternativeCall <- gsub("article46tvcXv2x_accountability", "article46tvc:v2x_accountability", model4AlternativeCall)
model4AlternativeCall <- as.formula(paste(model4AlternativeCall[2],
                                          model4AlternativeCall[1],
                                          model4AlternativeCall[3]))



model4AlternativeCall <- coxph(model4AlternativeCall, 
                               model4Data,
                               model4Data$genetic.any.weights)



model4interaction <- coxsimInteract(model4AlternativeCall,
                                    b1 = "article46tvc",
                                    b2 = "v2x_accountability",
                                    qi = "Marginal Effect", 
                                    X2 = seq(min(model4Data$v2x_accountability), 
                                             max(model4Data$v2x_accountability), 
                                             length.out = 100),
                                    means = FALSE, 
                                    expMarg = FALSE, 
                                    nsim = 1000, 
                                    ci = 0.95, 
                                    spin = FALSE,
                                    extremesDrop = TRUE)


model4interactionPlot <- simGG(model4interaction, 
                               xlab = "Government accountability", 
                               ylab = "Coefficient for remedial indications",
                               lcolour = "black", 
                               pcolour = "grey")+
  theme_classic()+
  geom_hline(aes(yintercept = 0), 
             linetype = "dashed", 
             colour = "grey")

model4interactionPlot$data$Lower50 <- model4interactionPlot$data$Min
model4interactionPlot$data$Upper50 <- model4interactionPlot$data$Max



pdf("model4interactionPlot.pdf")
print(model4interactionPlot)
dev.off()







##### Robustness checks: ####
robustness1 <- coxph(Surv(tstart, tstop, compliance) ~ 
                       article46tvc+
                       legislation +
                       jurisprudence +
                       exe_adm +
                       publication +
                       practical +
                       property +
                       reopen +
                       other_individual +
                       numberOfArticlesViolated+
                       years_since_legislative_election_tvc +
                       followCases+
                       polconiii+
                       nondemocracy+
                       newdemocracy1+
                       execap +
                       v2x_accountability+
                       year +
                       postProtocol11+
                       postCoMWorkingmethodChange +
                       postProtocol14 +
                       article2violated+
                       article3violated+
                       article5violated+
                       article6violated+
                       article8violated+
                       article10violated+
                       article13violated+
                       article14violated+
                       protocol1_1violated+
                       strata(respondent) +
                       cluster(respondent),
                     data = implementation,
                     x = TRUE)

summary(robustness1)

## This takes a while to run. Commented out to save time: ###
# genetic.any2 <- matchit(article46judgmentAnyPoint ~
#                           legislation +
#                           jurisprudence +
#                           exe_adm +
#                           publication +
#                           practical +
#                           prosecution +
#                           property +
#                           reopen +
#                           other_individual+
#                           numberOfArticlesViolated+
#                           article2violated+
#                           article3violated+
#                           article5violated+
#                           article6violated+
#                           article8violated+
#                           article10violated+
#                           article13violated+
#                           article14violated+
#                           protocol1_1violated+
#                           polconiii+
#                           nondemocracy+
#                           newdemocracy1+
#                           execap +
#                           v2x_accountability +
#                           Judgment.Year +
#                           postProtocol11+
#                           postCoMWorkingmethodChange +
#                           postProtocol14,
#                         data = matchingData,
#                         method = "genetic",
#                         unif.seed = 232321,
#                         wait.generations= 10,
#                         nboots = 1000,
#                         ratio = 2,
#                         pop.size = 1000,
#                         replace = TRUE,
#                         ties = TRUE)
# 
# save(genetic.any2,
#      file = "secondMatching.RData")
load("secondMatching.RData") ## loading results from second genetic matching




matchingData$genetic.any.weights2  <- genetic.any2$weights
implementation <- merge(implementation,
                        matchingData[,c("id",
                                        "genetic.any.weights2")],
                        by = "id",
                        all.x = TRUE)




robustness2 <- coxph(formula = Surv(tstart, tstop, compliance) ~ article46tvc + 
                       legislation +
                       jurisprudence +
                       exe_adm +
                       publication +
                       practical +
                       property +
                       reopen +
                       other_individual +
                       numberOfArticlesViolated+
                       years_since_legislative_election_tvc +
                       followCases+
                       polconiii+
                       nondemocracy+
                       newdemocracy1+
                       execap +
                       v2x_accountability +
                       cluster(respondent),
                     data = implementation[which(implementation$genetic.any.weights2 > 0), ], 
                     weights = implementation[which(implementation$genetic.any.weights2 > 0), ]$genetic.any.weights2,
                     y = TRUE,
                     x = TRUE,
                     tt = logtime)




matchingData2 <- na.omit(firstRow[,c('id',
                                     'article46judgmentAnyPoint',
                                     'legislation',
                                     'jurisprudence',
                                     'exe_adm',
                                     'publication',
                                     'practical',
                                     'prosecution',
                                     'property',
                                     'reopen',
                                     'other_individual',
                                     'numberOfArticlesViolated',
                                     "article2violated",
                                     "article3violated",
                                     "article5violated",
                                     "article6violated",
                                     "article8violated",
                                     "article10violated",
                                     "article13violated",
                                     "article14violated",
                                     "protocol1_1violated",
                                     "respondent",
                                     "Judgment.Year",
                                     "postProtocol11",
                                     "postCoMWorkingmethodChange",
                                     "postProtocol14")])

### Copied out to save time
# genetic.any3 <- matchit(article46judgmentAnyPoint ~
#                           legislation +
#                           jurisprudence +
#                           exe_adm +
#                           publication +
#                           practical +
#                           prosecution +
#                           property +
#                           reopen +
#                           other_individual+
#                           numberOfArticlesViolated+
#                           article2violated+
#                           article3violated+
#                           article5violated+
#                           article6violated+
#                           article8violated+
#                           article10violated+
#                           article13violated+
#                           article14violated+
#                           protocol1_1violated+
#                           Judgment.Year +
#                           postProtocol11+
#                           postCoMWorkingmethodChange +
#                           postProtocol14 +
#                           as.factor(respondent),
#                         data = matchingData2,
#                         exact =c(rep(F, times =23), rep(T, times = 47)),
#                         method = "genetic",
#                         unif.seed = 232321,
#                         wait.generations= 10,
#                         nboots = 1000,
#                         ratio = 1,
#                         pop.size = 1000,
#                         replace = TRUE,
#                         ties = TRUE)
# save(genetic.any3,
#      file = "thirdMatching.RData")

load("thirdMatching.RData")



summary(genetic.any3)

matchingData2$genetic.any.weights3 <- genetic.any3$weights





implementation <- merge(implementation,
                        matchingData2[,c("id","genetic.any.weights3")],
                        by = "id",
                        all.x = TRUE)

robustness3 <- coxph(formula = Surv(tstart, tstop, compliance) ~ article46tvc + 
                       legislation + jurisprudence + exe_adm + publication + practical + 
                       property + reopen + other_individual + numberOfArticlesViolated + 
                       followCases + 
                       strata(respondent) +
                       cluster(respondent), 
                     data = implementation[which(implementation$genetic.any.weights3 > 0), ], 
                     weights = implementation[which(implementation$genetic.any.weights3 > 0), ]$genetic.any.weights3, 
                     x = TRUE)

summary(robustness3)
cox.zph(robustness3, "rank")


robustness4 <- coxph(Surv(tstart, tstop, compliance) ~ 
                       article46tvc+
                       legislation +
                       jurisprudence +
                       exe_adm +
                       publication +
                       practical +
                       property +
                       reopen +
                       other_individual +
                       numberOfArticlesViolated+
                       years_since_legislative_election_tvc +
                       change_government_party+
                       followCases+
                       polconiii+
                       nondemocracy+
                       newdemocracy1+
                       execap +
                       v2x_accountability+
                       cluster(respondent),
                     weights = implementation[which(implementation$genetic.any.weights > 0),]$genetic.any.weights,
                     data = implementation[which(implementation$genetic.any.weights > 0),],
                     x = TRUE)

summary(robustness4)
cox.zph(robustness4, "rank")
robustness5 <- coxph(Surv(tstart, tstop, compliance) ~ 
                       article46tvc+
                       legislation +
                       jurisprudence +
                       exe_adm +
                       publication +
                       practical +
                       property +
                       reopen +
                       other_individual +
                       numberOfArticlesViolated+
                       years_since_legislative_election_tvc +
                       change_government_orientation +
                       followCases+
                       polconiii+
                       nondemocracy+
                       newdemocracy1+
                       execap +
                       v2x_accountability+
                       cluster(respondent),
                     weights = implementation[which(implementation$genetic.any.weights > 0),]$genetic.any.weights,
                     data = implementation[which(implementation$genetic.any.weights > 0),],
                     x = TRUE)

summary(robustness5)


robustnessChecks <- list(robustness1, 
                         robustness2,
                         robustness3, 
                         robustness4, 
                         robustness5)

for(i in 1:length(robustnessChecks)){
  robustnessChecks[[i]]$n <- sum(robustnessChecks[[i]]$y[,1] == 0)
}

texreg(robustnessChecks , 
       fontsize = "scriptsize",
       label = "tab:robustness",
       caption = "Robustness tests", 
       file = "robustnessChecks.tex",
       custom.model.names = paste("Model A", 1:length(robustnessChecks), sep = ""),
       include.rsquared = FALSE,
       include.maxrs= FALSE, 
       caption.above = TRUE, 
       include.missings = FALSE,
       include.zph = FALSE,
       custom.coef.map = list(article46tvc = "Remedial indications",
                              v2x_accountability  = "Government accountability", 
                              years_since_legislative_election_tvc  = "Years since last election",
                              change_government_party = "Change in government (party)",
                              change_government_orientation = "Change in government orientation",
                              polconiii =  "Political constraints",
                              nondemocracy = "Non-democracy",
                              newdemocracy1 = "New democracy",
                              execap = "Bureaucratic capacity",
                              legislation =  "Need for legislative change",
                              jurisprudence = "Need for jurisprudential change",
                              exe_adm = "Need for executive action",
                              publication = "Need for publication/dissemination",
                              practical = "Need for practical measures",
                              property =  "Need for property return",
                              reopen =  "Need for reopening of domestic case",
                              other_individual =  "Need for other individual measure",
                              numberOfArticlesViolated =  "Number of articles violated",
                              followCases = "Repetitive cases",
                              year = "Linear time trend", 
                              postProtocol11  = "After protocol 11", 
                              postCoMWorkingmethodChange  = "After change in CoM working methods",
                              postProtocol14 = "After protocol 14",
                              article2violated  = "Right to life violation",
                              article3violated  = "Prohibition of torture violation", 
                              article5violated  = "Right to liberty violation",
                              article6violated  = "Right to fair trial violation",  
                              article8violated  = "Right to privacy and family life violation", 
                              article10violated = "Freedom of expression violation",
                              article13violated =  "Right to effective remedy violation",
                              article14violated = "Prohibition of discrimination violation",
                              protocol1_1violated = "Property rights violation"))


#### Additional robustness tests using different subcomponents of government accountability ####

accountabilityModelsData <- implementation[which(implementation$genetic.any.weights > 0),]
model4vertical <- coxph(Surv(tstart, tstop, compliance) ~ 
                          article46tvc * v2x_veracc +
                          legislation +
                          jurisprudence +
                          exe_adm +
                          publication +
                          practical +
                          property +
                          reopen +
                          other_individual +
                          numberOfArticlesViolated+
                          years_since_legislative_election_tvc +
                          followCases+
                          polconiii+
                          nondemocracy+
                          newdemocracy1+
                          execap +
                          cluster(respondent),
                        weights = accountabilityModelsData $genetic.any.weights,
                        data = accountabilityModelsData ,
                        tt = logtime,
                        x = TRUE)
summary(model4vertical)


interactionVertical <- simGG(coxsimInteract(model4vertical,
                                            b1 = "article46tvc", 
                                            b2 = "v2x_veracc",
                                            X2 = seq(min(accountabilityModelsData$v2x_veracc), max(accountabilityModelsData $v2x_veracc), length.out = 100),
                                            expMarg = FALSE),
                             lcolour = "black",
                             pcolour = "grey",
                             xlab = "Vertical Accountability", 
                             ylab = "Coefficient for\nremedial indications")+
  theme_classic()+
  geom_hline(aes(yintercept = 0), colour = "grey")
interactionVertical$data$Lower50 <- interactionVertical$data$Min
interactionVertical$data$Upper50 <- interactionVertical$data$Max

model4horizontal <- coxph(Surv(tstart, tstop, compliance) ~ 
                            article46tvc * v2x_horacc +
                            legislation +
                            jurisprudence +
                            exe_adm +
                            publication +
                            practical +
                            property +
                            reopen +
                            other_individual +
                            numberOfArticlesViolated+
                            years_since_legislative_election_tvc +
                            followCases+
                            polconiii+
                            nondemocracy+
                            newdemocracy1+
                            execap +
                            cluster(respondent),
                          weights = accountabilityModelsData$genetic.any.weights,
                          data = accountabilityModelsData,
                          tt = logtime,
                          x = TRUE)

interactionHorizontal <- simGG(coxsimInteract(model4horizontal,
                                              b1 = "article46tvc", 
                                              b2 = "v2x_horacc",
                                              X2 = seq(min(accountabilityModelsData$v2x_horacc), max(accountabilityModelsData$v2x_horacc), length.out = 100),
                                              expMarg = FALSE),
                               lcolour = "black",
                               pcolour = "grey",
                               xlab = "Horizontal Accountability", 
                               ylab = "Coefficient for\nremedial indications")+
  theme_classic()+
  geom_hline(aes(yintercept = 0), colour = "grey")

interactionHorizontal$data$Lower50 <- interactionHorizontal$data$Min
interactionHorizontal$data$Upper50 <- interactionHorizontal$data$Max

model4diagonal <- coxph(Surv(tstart, tstop, compliance) ~ 
                          article46tvc * v2x_diagacc +
                          legislation +
                          jurisprudence +
                          exe_adm +
                          publication +
                          practical +
                          property +
                          reopen +
                          other_individual +
                          numberOfArticlesViolated+
                          years_since_legislative_election_tvc +
                          followCases+
                          polconiii+
                          nondemocracy+
                          newdemocracy1+
                          execap +
                          cluster(respondent),
                        weights = implementation[which(implementation$genetic.any.weights > 0),]$genetic.any.weights,
                        data = implementation[which(implementation$genetic.any.weights > 0),],
                        tt = logtime,
                        x = TRUE)

interactionDiagonal<- simGG(coxsimInteract(model4diagonal,
                                           b1 = "article46tvc", 
                                           b2 = "v2x_diagacc",
                                           X2 = seq(min(accountabilityModelsData$v2x_diagacc),
                                                    max(accountabilityModelsData$v2x_diagacc), length.out = 100),
                                           expMarg = FALSE),
                            lcolour = "black",
                            pcolour = "grey",
                            xlab = "Diagonal Accountability", 
                            ylab = "Coefficient for\nremedial indications")+
  theme_classic()+
  geom_hline(aes(yintercept = 0), colour = "grey")

interactionDiagonal$data$Lower50 <- interactionDiagonal$data$Min
interactionDiagonal$data$Upper50 <- interactionDiagonal$data$Max



pdf("AdditionalAccountabilityInteractionPlots.pdf")
print(ggarrange(interactionVertical, interactionHorizontal, interactionDiagonal, ncol = 1, nrow = 3))
dev.off()


#### Sensitivity tests for interaction model

model4Data$InTails <- ifelse(model4Data$v2x_accountability < 0.066 | model4Data$v2x_accountability > 1.779,1,0)
model4Data$InTails <- ave(model4Data$InTails, 
                          model4Data$id, 
                          FUN = function(x){ifelse(sum(x)>0,1,0)})
Model4DataWithoutTails <- filter(model4Data,
                                 InTails == 0)

model4withoutTails <- coxph(Surv(tstart, tstop, compliance) ~ 
                              article46tvc *  v2x_accountability +
                              legislation +
                              jurisprudence +
                              exe_adm +
                              publication +
                              practical +
                              property +
                              reopen +
                              other_individual +
                              numberOfArticlesViolated+
                              years_since_legislative_election_tvc +
                              followCases+
                              polconiii+
                              nondemocracy+
                              newdemocracy1+
                              execap +
                              cluster(respondent),
                            weights = Model4DataWithoutTails$genetic.any.weights,
                            data = Model4DataWithoutTails, 
                            tt = logtime,
                            x = TRUE)

summary(model4withoutTails)

model4interactionNoTails <- simGG(coxsimInteract(model4withoutTails,
                                                             b1 = "article46tvc",
                                                             b2 = "v2x_accountability",
                                                             qi = "Marginal Effect", 
                                                             X2 = seq(min(Model4DataWithoutTails$v2x_accountability), 
                                                                      max(Model4DataWithoutTails$v2x_accountability), 
                                                                      length.out = 100),
                                                             means = FALSE, 
                                                             expMarg = FALSE, 
                                                             nsim = 1000, 
                                                             ci = 0.95, 
                                                             spin = FALSE,
                                                             extremesDrop = TRUE),
                                              lcolour = "black", 
                                              pcolour = "grey")+
  theme_classic()+
  geom_hline(aes(yintercept = 0), linetype= "dashed")+
  xlab("Government Accountability") +
  ylab("Coefficient for Remedial Indications")

model4interactionNoTails$data$Lower50 <- model4interactionNoTails$data$Min
model4interactionNoTails$data$Upper50 <- model4interactionNoTails$data$Max

pdf("model4interactionNoTails.pdf")
print(model4interactionNoTails)
dev.off()

model4Data$influentialObs <- ifelse(abs(residuals(model4, "dfbetas")[,3]) > 2/sqrt(model4$n),1,0)
model4Data$influentialSpell <- ave(model4Data$influentialObs, 
                                   model4Data$id, 
                                   FUN = function(x){ifelse(sum(x)>1,1,0)})

Model4dataNoInfluentialSpells <- filter(model4Data, 
                                        influentialSpell == 0)

model4NoInfluentialSpells <- coxph(formula = Surv(tstart, tstop, compliance) ~ article46tvc * v2x_accountability  + legislation + 
                                     jurisprudence + exe_adm + publication + practical + property + 
                                     reopen + other_individual + numberOfArticlesViolated + years_since_legislative_election_tvc + 
                                     followCases + polconiii + nondemocracy + newdemocracy1 + 
                                     execap + cluster(respondent), data = Model4dataNoInfluentialSpells,
                                   weights = Model4dataNoInfluentialSpells$genetic.any.weights, 
                                   x = TRUE, tt = logtime)

summary(model4NoInfluentialSpells)
model4NoInfluentialSpells$n

model4interactionNoInfluentialSpells <- simGG(coxsimInteract(model4NoInfluentialSpells,
                                                             b1 = "article46tvc",
                                                             b2 = "v2x_accountability",
                                                             qi = "Marginal Effect", 
                                                             X2 = seq(min(Model4dataNoInfluentialSpells$v2x_accountability), 
                                                                      max(Model4dataNoInfluentialSpells$v2x_accountability), 
                                                                      length.out = 100),
                                                             means = FALSE, 
                                                             expMarg = FALSE, 
                                                             nsim = 1000, 
                                                             ci = 0.95, 
                                                             spin = FALSE,
                                                             extremesDrop = TRUE),
                                              lcolour = "black", 
                                              pcolour = "grey")+
  theme_classic()+
  geom_hline(aes(yintercept = 0), linetype= "dashed")+
  xlab("Government Accountability") + 
  ylab("Coefficient for Remedial Indications")

model4interactionNoInfluentialSpells$data$Lower50 <- model4interactionNoInfluentialSpells$data$Min
model4interactionNoInfluentialSpells$data$Upper50 <- model4interactionNoInfluentialSpells$data$Max

pdf("model4interactionNoInfluentialSpells.pdf")
print(model4interactionNoInfluentialSpells)
dev.off()



#### Tables for additional interaction models 

additionalInteractionModels <- list(model4vertical, model4horizontal, model4diagonal, model4withoutTails, model4NoInfluentialSpells)


for(i in 1:length(additionalInteractionModels)){
  additionalInteractionModels[[i]]$n <- sum(additionalInteractionModels[[i]]$y[,1] == 0)
}

texreg(additionalInteractionModels, 
       fontsize = "scriptsize",
       label = "tab:robustnessInteraction",
       caption = "Robustness tests for Model 4 of the letter", 
       file = "additionalInteractionModels.tex",
       custom.model.names = paste("Model A", 6:(length(robustnessChecks)+5), sep = ""),
       include.rsquared = FALSE,
       include.maxrs= FALSE, 
       caption.above = TRUE, 
       include.missings = FALSE,
       include.zph = FALSE,
        custom.coef.map = list(article46tvc = "Remedial indications",
                               `v2x_veracc` = "Vertical Accountability",
                               `article46tvc:v2x_veracc` = "Vertical Accountability * Remedial indications",
                               `v2x_horacc` = "Horizontal Accountability", 
                               `article46tvc:v2x_horacc` = "Horizontal Accountability * Remedial indications", 
                               v2x_diagacc = "Diagonal Accountability", 
                               `article46tvc:v2x_diagacc` = "Diagonal Accountability * Remedial indications",
                               v2x_accountability  = "Government accountability", 
                               article46tvcXv2x_accountability = "Government accountability * Remedial indications",
                               `article46tvc:v2x_accountability` = "Government accountability * Remedial indications",
                              years_since_legislative_election_tvc  = "Years since last election",
                              change_government_party = "Change in government (party)",
                              change_government_orientation = "Change in government orientation",
                              polconiii =  "Political constraints",
                              nondemocracy = "Non-democracy",
                              newdemocracy1 = "New democracy",
                              execap = "Bureaucratic capacity",
                              legislation =  "Need for legislative change",
                              jurisprudence = "Need for jurisprudential change",
                              exe_adm = "Need for executive action",
                              publication = "Need for publication/dissemination",
                              practical = "Need for practical measures",
                              property =  "Need for property return",
                              reopen =  "Need for reopening of domestic case",
                              other_individual =  "Need for other individual measure",
                              numberOfArticlesViolated =  "Number of articles violated",
                              followCases = "Repetitive cases",
                              year = "Linear time trend",
                              postProtocol11  = "After protocol 11",
                              postCoMWorkingmethodChange  = "After change in CoM working methods",
                              postProtocol14 = "After protocol 14",
                              article2violated  = "Right to life violation",
                              article3violated  = "Prohibition of torture violation",
                              article5violated  = "Right to liberty violation",
                              article6violated  = "Right to fair trial violation",
                              article8violated  = "Right to privacy and family life violation",
                              article10violated = "Freedom of expression violation",
                              article13violated =  "Right to effective remedy violation",
                              article14violated = "Prohibition of discrimination violation",
                              protocol1_1violated = "Property rights violation"))



###### Additional interaction effects for the appendix #####



interaction1 <- coxph(formula = Surv(tstart, tstop, compliance) ~ article46tvc * polconiii + 
                        v2x_accountability + legislation + 
                        jurisprudence + exe_adm + publication + practical + property + 
                        reopen + other_individual + numberOfArticlesViolated + years_since_legislative_election_tvc + 
                        followCases + polconiii + nondemocracy + newdemocracy1 + 
                        execap + cluster(respondent), data = model4Data, weights = model4Data$genetic.any.weights, 
                      x = TRUE)

interaction1plot <- simGG(coxsimInteract(interaction1,
                                         b1 = "article46tvc", 
                                         b2 = "polconiii",
                                         X2 = seq(min(model4Data$polconiii), max(model4Data$polconiii), length.out = 100),
                                         expMarg = FALSE),
                          lcolour = "black",
                          pcolour = "grey",
                          xlab = "Political constraints", 
                          ylab = "Coefficient for remedial indications")+
  theme_classic()+
  geom_hline(aes(yintercept = 0), colour = "grey")
interaction1plot$data$Lower50 <- interaction1plot$data$Min
interaction1plot$data$Upper50 <- interaction1plot$data$Max


pdf("interaction1plot.pdf")
print(interaction1plot)
dev.off()


interaction2 <- coxph(formula = Surv(tstart, tstop, compliance) ~ article46tvc * execap + 
                        v2x_accountability + legislation + 
                        jurisprudence + exe_adm + publication + practical + property + 
                        reopen + other_individual + numberOfArticlesViolated + years_since_legislative_election_tvc + 
                        followCases + polconiii + nondemocracy + newdemocracy1 + 
                        execap + cluster(respondent), data = model4Data, weights = model4Data$genetic.any.weights, 
                      x = TRUE, tt = logtime)

summary(interaction2)

interaction2plot <- simGG(coxsimInteract(interaction2,
                                         b1 = "article46tvc", 
                                         b2 = "execap",
                                         X2 = seq(min(model4Data$execap), max(model4Data$execap), length.out = 100),
                                         expMarg = FALSE),
                          lcolour = "black",
                          pcolour = "grey",
                          xlab = "Bureaucratic capacity", 
                          ylab = "Coefficient for remedial indications")+
  theme_classic()+
  geom_hline(aes(yintercept = 0), colour = "grey")
interaction2plot$data$Lower50 <- interaction2plot$data$Min
interaction2plot$data$Upper50 <- interaction2plot$data$Max

pdf("interaction2plot.pdf")
print(interaction2plot)
dev.off()


interactions <- list(interaction1, 
                     interaction2)

for(i in 1:length(interactions)){
  interactions[[i]]$n <- sum(interactions[[i]]$y[,1] == 0)
}

texreg(interactions , 
       fontsize = "scriptsize",
       label = "tab:interactions",
       caption = "Additional interaction models", 
       file = "interactions.tex",
       custom.model.names = paste("Model A", 11:12, sep = ""),
       include.rsquared = FALSE,
       include.maxrs= FALSE, 
       caption.above = TRUE, 
       include.missings = FALSE,
       include.zph = FALSE,
       custom.coef.map = list(article46tvc = "Remedial indications",
                              polconiii =  "Political constraints",
                              `article46tvc:polconiii` = "Remedial indications*Political constraints", 
                              execap = "Bureaucratic capacity",
                              `article46tvc:execap` = "Remedial indications*Bureaucratic capacity", 
                                v2x_accountability  = "Government accountability",
                              years_since_legislative_election_tvc  = "Years since last election",
                              change_government_party = "Change in government (party)",
                              change_government_orientation = "Change in government orientation",
                              nondemocracy = "Non-democracy",
                              newdemocracy1 = "New democracy",
                              legislation =  "Need for legislative change",
                              jurisprudence = "Need for jurisprudential change",
                              exe_adm = "Need for executive action",
                              publication = "Need for publication/dissemination",
                              practical = "Need for practical measures",
                              property =  "Need for property return",
                              reopen =  "Need for reopening of domestic case",
                              other_individual =  "Need for other individual measure",
                              numberOfArticlesViolated =  "Number of articles violated",
                              followCases = "Repetitive cases",
                              year = "Linear time trend", 
                              postProtocol11  = "After protocol 11", 
                              postCoMWorkingmethodChange  = "After change in CoM working methods",
                              postProtocol14 = "After protocol 14",
                              article2violated  = "Right to life violation",
                              article3violated  = "Prohibition of torture violation", 
                              article5violated  = "Right to liberty violation",
                              article6violated  = "Right to fair trial violation",  
                              article8violated  = "Right to privacy and family life violation", 
                              article10violated = "Freedom of expression violation",
                              article13violated =  "Right to effective remedy violation",
                              article14violated = "Prohibition of discrimination violation",
                              protocol1_1violated = "Property rights violation"))





### The End ####




